OpenHDF5
Kristin Braziunas
June 20, 2016
Learning to open HDF5 files in R
# load libraries
library(raster)## Loading required package: splibrary(rhdf5)
library(rgdal)## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3Let’s open a file
# set file path
f <- "../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"
# view h5 structure
h5ls(f)## group name otype dclass
## 0 / ATCOR_Input_File H5I_DATASET STRING
## 1 / ATCOR_Processing_Log H5I_DATASET STRING
## 2 / Aerosol Optical Depth H5I_DATASET INTEGER
## 3 / Aspect H5I_DATASET FLOAT
## 4 / Cast Shadow H5I_DATASET INTEGER
## 5 / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6 / Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7 / Illumination Factor H5I_DATASET INTEGER
## 8 / Path Length H5I_DATASET FLOAT
## 9 / Processing Version H5I_DATASET STRING
## 10 / Reflectance H5I_DATASET INTEGER
## 11 / Shadow_Processing_Log H5I_DATASET STRING
## 12 / Sky View Factor H5I_DATASET INTEGER
## 13 / Skyview_Processing_Log H5I_DATASET STRING
## 14 / Slope H5I_DATASET FLOAT
## 15 / Slope_Aspect_Processing_Log H5I_DATASET STRING
## 16 / Solar Azimuth Angle H5I_DATASET FLOAT
## 17 / Solar Zenith Angle H5I_DATASET FLOAT
## 18 / Surface Elevation H5I_DATASET FLOAT
## 19 / Visibility Index Map H5I_DATASET INTEGER
## 20 / Water Vapor Column H5I_DATASET INTEGER
## 21 / coordinate system string H5I_DATASET STRING
## 22 / flightAltitude H5I_DATASET FLOAT
## 23 / flightHeading H5I_DATASET FLOAT
## 24 / flightTime H5I_DATASET FLOAT
## 25 / fwhm H5I_DATASET FLOAT
## 26 / map info H5I_DATASET STRING
## 27 / to-sensor azimuth angle H5I_DATASET FLOAT
## 28 / to-sensor zenith angle H5I_DATASET FLOAT
## 29 / wavelength H5I_DATASET FLOAT
## dim
## 0 1
## 1 1
## 2 544 x 578 x 1
## 3 544 x 578 x 1
## 4 544 x 578 x 1
## 5 544 x 578 x 1
## 6 544 x 578 x 1
## 7 544 x 578 x 1
## 8 544 x 578 x 1
## 9 1
## 10 544 x 578 x 426
## 11 1
## 12 544 x 578 x 1
## 13 1
## 14 544 x 578 x 1
## 15 1
## 16 1 x 1
## 17 1 x 1
## 18 544 x 578 x 1
## 19 544 x 578 x 1
## 20 544 x 578 x 1
## 21 1
## 22 5732053
## 23 5732053
## 24 5732053
## 25 426 x 1
## 26 1
## 27 544 x 578 x 1
## 28 544 x 578 x 1
## 29 426 x 1Import spatial information
# import spatial info
mapInfo <- h5read(f,
"map info",
read.attributes = TRUE)
mapInfo## [1] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## attr(,"Description")
## [1] "Basic Map information for envi style programs"Grab reflectance metadata
# read in reflectance data attributes
reflInfo <- h5readAttributes(f, "Reflectance")
reflInfo## $DIMENSION_LABELS
## [1] "Wavelength" "Line" "Sample"
##
## $Description
## [1] "Atmospherically corrected reflectance."
##
## $`Scale Factor`
## [1] 10000
##
## $Unit
## [1] "unitless. Valid range 0-1."
##
## $`data ignore value`
## [1] "15000.0"# define scale factor
scaleFactor <- reflInfo$`Scale Factor`
# define no data value
noDataValue <- reflInfo$`data ignore value` # it's a character
noDataValue <- as.numeric(noDataValue)
str(noDataValue) # looks good## num 15000Import data dimensions
# open file for viewing
fid <- H5Fopen(f)
# open the reflectance dataset
did <- H5Dopen(fid, "Reflectance")
did # note dimensions are in columns x rows x bands## HDF5 DATASET
## name /Reflectance
## filename
## type H5T_STD_I16LE
## rank 3
## size 544 x 578 x 426
## maxsize 544 x 578 x 426# grab dataset dimensions
sid <- H5Dget_space(did)
dims <- H5Sget_simple_extent_dims(sid)$size
dims # columns, rows, bands; r reads rows 1st, columns 2nd## [1] 544 578 426# close all open connections
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)Read in reflectance data
# extract slice of H5 file
b56 <- h5read(f,
"Reflectance",
index = list(1:dims[1], 1:dims[2], 56))
b56
class(b56)Convert data to matrix
# convert to matrix
b56 <- b56[,,1] # z dimension goes away, set to 1
class(b56)## [1] "matrix"# let's plot some data
image(b56)
# stretch the image by log transforming
image(log(b56), main = "Log transformed data")
# force non-scientific notation
options("scipen"=100, "digits"=4)
# look at histograms
hist(b56)
# histogram of stretched image
hist(log(b56))
Time to clean up our data
Note that data is stored as integers rather than floats. This way, data takes up less space. Can use scale factor to convert to decimal point at later time.
# assign no data values to object
b56[b56 == noDataValue] <- NA
# apply scale factor
b56 <- b56 / scaleFactor # grabbed scale factor earlier
hist(b56) # right skew
Transpose data
# transpose data, flip rows and columns
b56 <- t(b56)
image(log(b56)) # still not quite right
Create spatial extent
# split out map info object
mapInfo <- strsplit(mapInfo, ",")
mapInfo <- unlist(mapInfo) # remove nesting
mapInfo## [1] "UTM" "1" "1"
## [4] "325963.0" "4103482.0" "1.0000000000e+000"
## [7] "1.0000000000e+000" "11" "North"
## [10] "WGS-84" "units=Meters"# value at 3 element in list
mapInfo[3] # note that it's a character## [1] "1"# define upper left hand corner coordinate
xMin <- as.numeric(mapInfo[4])
yMax <- as.numeric(mapInfo[5])
# get spatial resolution
xRes <- as.numeric(mapInfo[6])
yRes <- as.numeric(mapInfo[7])
# define lower right hand corner coordinate
xMax <- xMin + (dims[1] * xRes)
yMin <- yMax - (dims[2] * yRes)
# create extent object
rasExt <- extent(xMin, xMax, yMin, yMax)
rasExt## class : Extent
## xmin : 325963
## xmax : 326507
## ymin : 4102904
## ymax : 4103482# create raster object
b56r <- raster(b56,
crs = CRS("+init=epsg:32611"))
extent(b56r) <- rasExt
b56r## class : RasterLayer
## dimensions : 578, 544, 314432 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4102904, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0.0033, 0.5678 (min, max)# plot data
plot(b56r, main="Spatially referenced data")
Import NEON functions
# install devtools
# install.packages("devtools")
library(devtools)
# install_github("lwasser/neon-aop-package/neonAOP")
library(neonAOP)
b55 <- open_band(f,
bandNum = 55,
epsg = 32611)
plot(b55)
# import several bands
bands <- c(58, 34, 19) # decreasing numeric order intentional to make r, g, b
bands## [1] 58 34 19# create raster stack
RGBStack <- create_stack(f,
bands = bands,
epsg = 32611)
plot(RGBStack) # plot 3 bands separately
# plot RGB image
plotRGB(RGBStack,
stretch = "lin") # make sure ordered correctly
# cir image
bands <- c(90, 34, 19)
# create raster stack
CIRStack <- create_stack(f,
bands = bands,
epsg = 32611)
# plot RGB image
plotRGB(CIRStack,
stretch = "lin") # make sure ordered correctly